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I ' Abstract 

^ : 

^ , We propose a computational method to obtain series expansions 

in powers of time for general dynamical systems described by a set 
of hierarchical rate equations. The method is generally applicable to 
• i-H . problems in both equilibrium and nonequilibrium statistical mechanics 

^ I such as random sequential adsorption, diffusion-reaction dynamics, 

■ and Ising dynamics. New result of random sequential adsorption of 

dimers on a square lattice is presented. 

PACS numbers: 05.50.+q,05.20.-y,05.70.Ln 

Most dynamical models in equilibrium and nonequilibrium statistical me- 
chanics can be described by a set of hierarchical rate equations: random 
sequential adsorptions (RSA) and their variants [IH, diffusion- reaction mod- 
els 0, and kinetic Ising models 0]. The specifications of interactions between 
the components of the system or interactions between the environment and 
the system give a deterministic time evolution for the distribution functions 
once the initial condition of the system is given. Exact solutions are often 
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restricted to simple models, while one has to use approximate methods for 
those that resist exact analyses. 

Power series expansion is one of the methods of controlled approxima- 
tions. There are a number of systematic methods to obtain time-power se- 
ries [§, 1^, 1^, 0]. Long series were usually obtained by reducing a problem 
to an enumeration problem similar to the series expansions in equilibrium 
statistical mechanics [§. 

Lattice enumeration problems are typically limited by CPU time due to 
their exponential growth. In this Letter, we discuss a general method to ob- 
tain series expansions based on rate equations for lattice models. Our method 
differs from previous approaches in that the limiting factor is memory space, 
but it is faster in time. The method is applicable to many different prob- 
lems, in particular RSA, RSA with diffusional relaxation, reaction-diffusion 
problems, and general Ising dynamics. 

To illustrate the method, we take the RSA of random dimer filling on a 
square lattice as an example P|. Dimers of random orientations are dropped 
randomly and sequentially, at a rate of k per lattice site per unit time, onto 
an initially empty, infinite square lattice. Hereafter we set k equal to unity 
without loss of generality. If the chosen two nearest neighbor sites are un- 
occupied, the dimer is adsorbed on the lattice. If one of the chosen sites is 
occupied by a previously adsorbed dimer, the adsorption attempt is rejected. 
The first few of an infinite number of rate equations for this process are given 

4P(oo), (1) 

P(oo) -2P(ooo) -4P(8o), (2) 
2P(ooo) - 2P(oooo) - 4P(8oo) - 2P(o8o), (3) 

2P(§o) - 2P(8§) - 2P(o8o) - 2P(8oo) - 2P(o8°), (4) 

where P(C) stands for the probability of finding a configuration C of sites 
specified empty or filled. An o denotes an empty site configuration. We 
have taken into account the symmetries of the problem, i.e. the invariance 
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of P{C) under all lattice group operations (e.g. rotation and reflection). For 
notational convenience we will use Pc = P{C). 

Let Cq denotes a particular configuration of interest, then Pc^ is the 
configuration probability associated with Cq. On physical grounds, we expect 
Pco to be a well behaved function of time t, and one would expect to obtain 

the Taylor expansion Pcoit) = J2 Pco t"'/^^-: with n-th derivative of Pcq 

n=0 



given by 



p (n) d'^Pco (t) 
^Co = 



dt"^ 



(5) 

t=o 



The zeroth derivative of P^o is determined by the initial condition. The 
first derivative of Pcq is obtained by the rate equation associated with Cq. 
For the random dimer deposition problem, we choose the configuration of 
an empty site to be Cq. We have, from Eq. (|I|) to (^), P(o)(°) = 1 and 
P(o)(i) = -4, P(oo)(i) = -7, P(ooo)« = -10, and P(8o)^^^ = -10. To 
compute the second derivative of P(o), we take the first derivative of Eq. ([^) 
and use the result of P(oo)(i) to obtain P(o)(2) = (_4)(_7) = 28. For 
the third derivative P(o)(^), we take the second derivative of Eq. (0), which 
in turn needs first derivative of Eq. (^. Using Eq. (2) to (4) for the first 
derivatives, we get P(o)(3) = —268. 

To computerize the calculations with high efficiency, we make the follow- 
ing important observations. For any configuration C, the rate equation is 
always of the form 

J2Xc'Pc'. (6) 



dPc 



dt ^, 



This immediately gives us the first derivative of Pc', in particular the first 
derivative of Pcq of interest. Since the derivative operator is linear, the second 
derivative of Pc is a linear combination of the first derivatives of Pc on the 
right-hand side of the rate equation @. For each higher derivative, new rate 
equations and new configurations are involved. Let Gi denote the set of new 
configurations generated in the calculation of the i-th derivative of Pcp, and 
Gl the corresponding j-th derivatives of the set of configurations. We observe 
that . . ., (determined at the {n - l)-th derivative), G^'^ 

G""^, . . ., G°_2 (determined at the (n — 2)-th derivative), . . ., Gq are known 
before calculating the n-th derivative. In other words, Gj are predetermined 
where 0<i + j<n — lat this stage. The derivatives Gq, G"~^, . . ., Gn_i, 
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would then be determined in a bottom-up fashion by recursive use of the 
rate equations. The method is efficient in the sense that each value in 
< i < n, is calculated exactly once and the rate equation for a configuration 
C is also generated exactly once. By repeated use of this procedure, we can 
in principle calculate Gq for any n > 1. 

In most cases, the knowledge of the initial condition of the system (usually 
one starts with an initially empty lattice) can further improve the efficiency 
of the method. For example, let h be the highest order of derivative of 
Pcf) for which we wish to compute. Then only a subset G/, i < h, of 
Gi needs to be considered, once it is known that an element in Gi — G/ 
gives a zero contribution in the subsequent calculations. This observation is 
particularly useful for models which introduce particles in the rate equations; 
this is typical in models involving desorption or diffusion. 

To implement the above scheme, we store symbolically the rate equations 
and the intermediate derivatives G^. For each configuration C, we store the 
configuration and its first j derivatives known so far. The rate equation for 
C is generated when the first derivative of Pc is needed, and the equation 
is represented as an array of pointers to other configurations corresponding 
to the right-hand side of Eq. (^) , and the corresponding array of coefficients 
Ac". 

The structural part implementing the recursive use of rate equations is 
independent of the details of a model. The function calls for a core routine to 
generate a rate equation if it is not available. Symmetries of the problem are 
treated here. This routine is called only once for each configuration C. To 
facilitate this, only a pointer of reference is used for each unique configuration. 
Configurations are searched with the help of a hash table or a tree data 
structure. 

The merit of the approach allows us to obtain a series up to for the 
RSA of dimers on a square lattice (Table |I]). Results for other models will 
be published elsewhere. The computational costs are presented in Table 
2. Poland obtained the same series up to [0]. For the RSA of monomer 
filling with nearest-neighbor exclusion on a square lattice, the efficiency of our 
algorithm is comparable to the algorithm based on the operator formalism 



[0], but it is not as good as the method used by Baram and Fixman pi 
However, the power of our method is its generality. 

The approach to the jamming state for lattice RSA is often exponential. 
We follow the analysis of series by Dickman, Wang, and Jensen Wj. First 
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we transform the coverage 6{t) = 1 — P{o^t) into a function of y = 1 — e~*. 
Another transformation via z = {1 — e^^^)/6, is performed for the second 
time, and we examine various Fade approximants to the z series. As we see 
from Fig. |l], the resulting estimates for 6^ are excellent for 1.736 < b < 1.742. 
The intersections of [6, 8], [7, 7], [7, 6], [8, 6], [6, 7], and [5, 9] approximants 
around this range of b yield an estimate of 9oo = 0.9068088(4), where the 
last digit denotes the uncertainty. This result is in good agreement with the 



simulation resuh of 0.906873 ± 0.000138 [g or the result of 0.9068 obtained 
via approximate truncation procedures [^. 

In summary, we introduced a new approach to deriving power series ex- 
pansions in time, and we have applied the method to a RSA problem. For 
the RSA of dimers on a square lattice, we are able to generate a 15-term 
series, thence to derive the most precise estimate for the jamming coverage 
yet presented. Our computational method is general and it can be used to 
handle a variety of problems based on rate equations, especially those deal 
with the kinetics of the lattice models. 

This work is supported in part by a National University of Singapore 
Academic Research Grant RF950601. Calculations were performed on the 
facilities of the Computation Center of the Institute of Fhysical and Chemical 
Research. 



5 



References 

[I] Evans J W 1993 Rev. Mod. Phys. 65 1281 

[2] Kuzovkov V and Kotomin E 1988 Rep. Prog. Phys. 51 1479 
[3] Glauber R J 1963 J. Math. Phys. 4 294 
[4] Evans J W 1984 Physica A 123 297 

[5] Poland D 1989 J. Stat. Phys. 55 389; 1990 J. Stat. Phys. 59 935 

[6] Baram A and Kutasov D 1989 J. Phys. A: Math. Gen. 22 L251 

[7] Dickman R, Wang J -S and Jensen I 1991 J. Chem. Phys. 94 8252 

[8] Martin J L 1974 Phase transitions and Critical Phenomena., Vol 3 eds 
C Domb and M S Green (New York: Academic Press) p97. 

[9] Nord R S and Evans J W 1985 J. Chem. Phys. 82 2795 

[10] Poland D 1991 Langmuir 7 514 

[II] Baram A and Fixman M 1995 J. Chem. Phys. 103 1929 
[12] Nord R S 1991 J. Statist. Comput. Simul. 39 231 



6 



Table 1: Taylor expansion coefficients for the probability of finding an empty 
site of random dimer filling on a square lattice. 



n 


P(o)(") 





1 


1 


-4 


2 


28 


3 


-268 


4 


3212 


5 


-45868 


6 


756364 


7 


-14094572 


8 


292140492 


9 


-6653993260 


10 


164952149516 


11 


-4416119044972 


12 


126863203272268 


13 


-3889473277203116 


14 


126677386324657804 
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Table 2: CPU time and memory usage on a DEC3000/900 for the dimer 
problem. 



order 


memory (megabytes) 


time (seconds) 


9 


1 


1.0 


10 


3 


4.4 


11 


12 


22.0 


12 


45 


114.7 


13 


173 


612.8 


14 


666 


3131.3 
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Figure 1: Fade approximant estimates for the jamming coverage ^oo as a 
function of the transformation parameter b, in the crossing region. 
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